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Polymer rheology simulations at the meso- and macroscopic scale 
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Abstract. - We show that simulations of polymer rheology at a fluctuating mesoscopic scale and 
at the macroscopic scale where flow instabilities occur can be achieved at the same time with 
dissipative particle dynamics (DPD) technique. We model the visco-elasticity of polymer liquids 
by introducing a flnite fraction of dumbbells in the standard DPD fluid. The stretching and 
tumbling statistics of these dumbbells is in agreement with what is known for isolated polymers in 
shear flows. At the same time, the model exhibits behaviour reminiscent of drag reduction in the 
turbulent state: as the polymer fraction increases, the onset of turbulence in plane Couette flow 
is pushed to higher Reynolds numbers. The method opens up the possibility to model nontrivial 
rheological conditions with ensuing coarse grained polymer statistics. 
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Introduction. — The plethora of intriguing phenom- 
ena that can be observed in flows of complex fluids is at- 
tracting increasing attention among physicists. The study 
of polymer fluids and melts has since long occupied a cen- 
tral position within this field PQH]. Typically, one is ei- 
ther interested in the response properties of polymeric flu- 
ids or the way they flow. The former set of problems 
concerns viscoelasticity of polymers and biomaterials, be- 
haviour of single long flexible molecules in flows, etcetera 
[5H1] . The latter usually focuses on the differences between 
ijiacroscopic flows of Newtonian and polymeric fluids in 
the same geometry. One such striking example is the re- 
cently discovered chaotic flows of dilute polymer solutions 
at very small Reynolds numbers - the so-called purely elas- 
tic turbulence 0[S]. Another is the phenomenon of drag- 
reduction - the observation that even minute amounts of 
polymer can signiflcantly suppress Newtonian turbulence 
and hence reduce turbulent drag [7] . In this Letter we de- 
velop a mesoscopic simulation method that is capable of 
addressing both classes of problems. 

Simulation methods of polymers essentially fall into two 
classes. On one hand there are many mesoscopic coarse- 
grained approaches that have mainly been developed to 
study the thermo- and hydrodynamic properties of poly- 
mers in equilibrium and in weakly non-equilibrium situa- 



tions such as an imposed small shear. Such results can, 
e.g., be compared with recent experimental results for the 
orientation statistics of single DNA molecules in solution 
[5]. However, these models typically cannot be scaled up 
to simulate flow at macroscopic rheological scales. On 
the other hand, numerical studies of polymer rheology at 
macroscopically relevant scales are based almost exclu- 
sively on numerical implementations of continuum con- 
stitutive equations like the Oldroyd-B or FENE-P models 
[T][5]. By their very nature, these deterministic approaches 
only give the average behaviour, so they cannot give in- 
sight into the coupling between the macroscopic flow be- 
haviour and the molecular properties. In this Letter, we 
for the first time bridge the gap between these two ap- 
proaches and scales by introducing a coarse-grained Dis- 
sipative Particle Dynamics (DPD) model [5HTU] for visco- 
elastic fiows which is a solution of elastic springs (dumb- 
bells). Unlike the previous DPD simulations of polymeric 
fluids jllj , we do not attempt to resolve internal dynamics 
of long polymer molecules. Instead we view the dumbbells 
as collective elastic degrees of freedom (normal modes) 
that render the solution viscoelastic. We show that the 
model is powerful enough to exhibit both stretching and 
tumbling of dumbbells at mesoscopic scales, and proper 
flow behaviour at hydrodynamic scales like the dramatic 
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polymer drag reduction of turbulence. 

Our approach holds an additional promise. Computa- 
tional rheology still turns out to be a major challenge [T^- 
The main difficulty lies in the fact that the convective non- 
linear term in the constitutive equation can lead to a lo- 
cal exponential growth of the components of the polymer 
stress tensor [13]. As a result, it is difficult to come up 
with robust numerical schemes. Although progress can 
be made by going to special variables [13], there is still 
a great need for an easy-to-use and robust method for 
simulating polymer rheology in complex geometries. As is 
well known, mesoscale models like DPD and Lattice Boltz- 
mann |14| are versatile methods to simulate laminar and 
turbulent Newtonian flows in rheometric flows |15] and in 
complicated geometries |16| : we expect our extension of 
the DPD model to have the same advantage. 

DPD model for a visco-elastic fluid. DPD [^[3] 

is an off-lattice method in which one simulates the dy- 
namics of particles which we can intuitively think of as 
representing mesoscopic blobs of fluid. The interparticle 
interactions consist of conserved, dissipative and random 
forces tuned to reproduce hydrodynamic behaviour on the 
scale of a few particles. Particles are assumed to have 
mass m and velocity Vi — Ri at positions Ri so that 
their equations of motion are 
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For the first three terms we follow [T5] and take the stan- 
dard choices [SlIS]: the force F'^°^^ is a soft repulsion. 
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Here Ri 



\Ri — Rj\ is the distance between the parti- 



cles, Rij the interparticle unit vector, and the coefficient a 
measures the strength of this interaction. The conserved 
force F™"*^ increases linearly when particles come within 
the range Rq. From here on we take m = Rq = 1, so that 
distances are in units of the particle radius. The dissipa- 
tive force F'^^^^ acts to equalise velocities of neighbouring 
particles, while the random force ir"^""^ represents a cou- 
pling to a heat bath. We use the standard form [SJin] 

Ff' = -T{R,, ■ V,,) max(l - i?,, , O)^/^^,;, , (3) 
Flf'" = ^2TT/M max( 1 - i?,, , 0) ^/^R,j , (4) 

where F is a friction coefficient, T is the temperature, 
Vij = Vi — Vj is the relative velocity vector. At is the 
integration timestep, and are independent Gaussian 
random variables with zero mean and unit variance. The 
form of iT'"''*"'^ is chosen to ensure that the fluctuation- 
dissipation theorem holds in the absence of applied flow 
PlfTT] . Note that the best convergence of the DPD model 
to hydrodynamic behaviour is achieved when all terms in 
Eq.(IT]) are roughly of the same order. 



model [TS] . It is motivated by the observation that the so- 
called Oldroyd-B constitutive equation is exact for non- 
interacting dumbbells with linear elastic springs [HE]- 
Linear Hookean springs, however, can lead to diverging 
stresses since the two partners can separate inflnitely far. 
We therefore build our DPD model for polymer flows by 
using nonlinear dumbbells: we assign a unique partner 
to a flnitc fraction of the particles and introduce, in the 
spirit of the Finitely Extensible Nonlinear Elastic (FENE) 
spring constitutive equations [Tl[2] , a FENE-force between 
each particle and its partner. 



■iFENE 



HfeneRi 



1 - (i?.,/i?n,ax)2 



(5) 



The FENE spring behaves as a Hookean spring with stiff- 
ness -ffpENE at small extensions and as a stiff rod when 
the extension Rij is close to its maximum value i?max- 

While this force works well in smooth flows, large lo- 
cal stresses in turbulent flows can, because of the finite 
timestep, sporadically lead to a blow-up due to configura- 
tions with Rij > i?max- It is therefore numerically advan- 
tageous to use instead in such situations a FENE-inspired 
nonlinear force that has no sharp maximum extension. 
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To the lowest order in the extension, this nonlinear force is 
the same as the FENE force ^ if we identify i?niax with 
i?c, but unlike Eq.© it remains finite for all Rij. We 
will use this soft force in the last section of this Letter. 
Once again, we stress that the dumbbells with the force 
law ([5]) or (|6]) do not represent individual polymers but 
rather collective viscoelastic degrees of freedom. 

The equations of motion are solved with a version of 
the velocity- Verlet algorithm [TH] . Unless noted otherwise, 
we take a = 1,T = 0.1, F = 2, particle density p = 4 
and timestep A< = 0.05. Simulations are performed in 
the plane Couette flow geometry with periodic boundary 
conditions in the x (streamwise) and y (spanwise) direc- 
tion, and no-slip walls perpendicular to the z (gradient) 
axis. The typical dimensions of our simulation box are 
LxXLyXLz = 20x20x20 so that when all particles are 
dumbbells we have 2 • (20)^ 16000 dumbbells at den- 
sity p — 4:. The walls are implemented as a soft repulsion 
potential in the normal direction — see jl8| for details. 

Rheological properties of the dumbbeU DPD 
modeL — Wc now show that this dumbbell DPD fluid 
exhibits the main characteristics of polymer rheology. In 
a simple shear flow, the stress tensor for an incompress- 
ible polymer solution <t has 3 independent components: 
there is one independent off-diagonal component while the 
diagonal components are linear combinations of the so- 
called flrst and second normal stress differences A^i and 
7V2- While these both vanish for a Newtonian fluid, the 
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normal stress differences characterise the rheological prop- 
erties of a viscoelastic fluid, in particular its relaxation 
time r''°'. More precisely, for our plane Couette geometry 
we have to lowest order in 7 = dv^/ dz, 

(Jxz = CFzx = -nil = \axx-'7zz \ = 277pt' 



(7) 



where 77 = rjs + "qp, with 775 and rjp being the solvent 
and the polymer contributions to the total viscosity. The 
expression for the first normal stress difference defines the 
effective polymer relaxation time t'"^^. The second normal 
stress difference N2 = \fJyy—<yzz\ Ni for most polymer 
solutions. For the Oldroyd-B constitutive equation, N2 
vanishes identically [Il[2]- 

To estimate the shear viscosity 77 and the relaxation time 
r"^' for different dumbbell parameters, we have performed 
simulations at various shear rates 7 and computed the 
components of the total stress tensor from the pairwise 
interactions by time-averaging the virial formula [TllISQ] : 



f 



LxLyL^ 
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over a time interval of length 200, which was long enough 
to ensure convergence; here 5Vl = Vi — Vj and is the 
average velocity of particles at the vertical position Zi. We 
have verified that only one shear stress component takes 



significant values (typically crxy/<7xz — o'yz/o'x 
and that axx ^ f — 



0.01), 

so that indeed N2 ^ A^i, as 
desired. The accuracy of ([5]) has also been validated by 
comparing its results to the total momentum transferred 
through the walls [18] . The viscosity 77 and the relaxation 
time r''*'' are then obtained from linear and quadratic fits 
of axz and A^i to 7 as illustrated in Fig. [TJa). By vary- 
ing the fraction fp of dumbbells from 0% to 100%, we 
can separate the individual contributions of rjs and ?7p; as 
shown in Fig.[TJb), the total viscosity increases essentially 
linearly with the concentration, so that for the choosen 
DPD parameters 775 « 0.68 and 77P ~ 0.29/p. 

At fixed polymer concentration, r''°' is approximately 
inversely proportional to -ffpENE for large -Rmax, as can be 
inferred from Table [T] This is as expected for the relax- 
ation of a dumbbell in a viscous medium [1]. (We have 
found that we can increase r''^' while keeping -ffpENE — 1 
by using chains with more than two beads [20] ) . We have 
checked that there is little shear-thinning in our model: 
T] decreases by at most 10% at high shear rates [l^; we 
have also observed that the magnitude of shear-thinning 
increases with t''°'. We have also estimated r'''^' from the 
decay of the shear stress autocorrelation function [TT] . As 
shown in Fig. [TJc), it can be approximated by a single 
exponential that yields r'''^' = 1.51 and t'"^^ = 6.1 for 
^max = 4, and iJpENE = 0.5 and i/pENE = 0.1, respec- 
tively. These results are within 20% of the values quoted 
in Tabled! 

The above results show that the dumbbell DPD model 
is a faithful mesoscopic representation of the Oldroyd-B 
type constitutive equation for polymer rheology. 
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0.77 4.18 
0.67 3.07 
0.49 2.01 
0.27 1.2 
0.15 0.64 


0.91 5.29 
0.76 3.55 
0.55 2.25 
0.29 1.23 
0.16 0.66 


1.04 5.63 
0.82 4.22 
0.58 2.72 
0.32 1.23 
0.16 0.65 



Table 1: Zero-shear rate fluid viscosity and relaxation time for 
various dumbbell parameters with fp — 100% and the dumb- 
bell FENE force ©. 



Polymer stretching and orientation statistics. 

The statistics of single-polymer stretching and orientation 
in shear flows has recently been studied both experimen- 
tally [3] and theoretically [21]. We now demonstrate that 
our DPD dumbbell model is also capable of qualitatively 
reproducing the single- molecule results {511^. 

We first study the extension of the FENE dumbbells 
(0 for various values of the Weissenberg number Wi = 
7T''''' that describes the balance between the shear flow 
and elasticity. Fig.|2ja) shows the Probability Distribution 
Function (PDF) of the dumbbell extension for three values 
of Wi. Upon the increase of the Weissenberg number, 
the shape of P{R/Rma,x) changes very much in line with 
what has been found for single polymers in shear flows by 
other methods [4l|2Tl|22]: a peak at small and a power- 
law decay at large extensions for Wi < 1, and a rather 
flat distribution across a large range of extensions above 
the coil-stretch transition for a single polymer at Wi « 1. 
At very high Wi the distribution becomes strongly peaked 
close to i?max but we have not explored this regime. 

The average orientation of a polymer in a shear flow, 
which in our dumbbell model translates into the average 
orientation of the dumbbell, is characterised by the spher- 
ical angles d and </>. The distribution of 9, the angle be- 
tween the dumbbell and the xz shear plane, is found to 
be well approximated by a Lorentzian, in agreement with 
single-dumbbell models [2^11^ . We focus here on the PDF 
of the angle (j) that the projection of each dumbbell in the 
shear plane makes with the flow direction; in view of the 
symmetry of the system, we only need to consider angles 
-7r/2 < 4>< tt/2. 

The distribution P(0) is shown in Fig. ^h) together 
with the flts to the expression 



P(0) C 
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-(/? + sin (f>j dip, 
(9) 

which can be derived from the Fokkcr-Planck equation for 
P{(j),t) [22] • In this expression D plays the role of a dif- 
fusion coefficient for the orientation angle (j). As can be 
seen from the figure, the above expression fits our data 
very well, and this allows us to study the shear-rate de- 
pendence of D; from the fitted values of j/D listed in the 
caption, one infers that D decreases with increasing shear 
rate approximately as 7"^ '^. Although for a FENE model 
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Fig. 1: (colour online) (a) The normal stress difference A''i as a function of 7^ for our DPD model ^ with ffpENE = 0.5 and 
i?max = 4 and various dumbbell fractions fp. r"^' is determined from the fit of A'^i to 7^. (b) Viscosity ratio r]p/i]s as a function 
of polymer fraction fp for two sets of dumbbell parameters, (c) The normalized (dumbbells contribution to the) shear stress 
autocorrelation function as a function of time. Lines are single-exponential fits: red - 1.31 e"*''^'^'^, blue - 1.0 e~ 



-t/6.1 










i/pENE 
^tum 


0.5 1 

29.80 20.36 


0.5 1 

29.20 23.01 



Fig. 2: (colour online) (a) Extension PDF for Wi — 0.5,2 
and 5 and Hfene ~ 0.5 and i?max ~ 5.0. For Wi > 1 the 
tail of the PDF is no longer algebraic [21]; this is a way to 
characterise the coil-stretch transition (Wi « 1). (b) PDF of 
(j> for a DPD polymer with J/fene = 1-0 and i?max ~ 5.0 for 
different Weissenberg numbers. Continuous lines are fits to 
From low to high Wi, yellow: j/D — 2.9, red: j/D — 41.6, 
blue: j/D = 97.6, orange: 7/D = 322.9. 



one does not expect a pure scaling form, we can compare 
this to what is expected for a single dumbbell with a linear 
spring, for which D ^ 7^^, as well as to what is expected 
for inextensible rods, for which D = const [H]. This com- 
parison indicates that over the range of shear rates con- 
sidered, the finite extensibility of the dumbbells already 
reduces the shear rate dependence of D noticeably. 

Tumbling statistics. — As we noted before, the ori- 
entation statistics of an isolated polymer can largely be 
understood by studying a dumbbell in a shear flow [3T]. 
This also holds for the tumbling statistics [511^ : in a shear 



Table 2: Single dumbbell tumbling times at Wi = 1 for various 
dumbbells stiffness and maximal extension. Essentially the 
same values are found for all Weissenberg numbers up to Wi ~ 
5, except that we find a slight decrease with increasing shear 
rate for the parameters corresponding to the rightmost column. 



flow, the polymer is most of the time almost aligned with 
the mean flow, but due to the torque exerted on it by 
shear, every now and then it rapidly tumbles before com- 
ing back close to the average orientation; the tumbling 
time T*"™ is defined as the typical time between two such 
events [511^ . In our DPD model we have many dumbbells 
which interact via Nevertheless, we shall now show 

that our model still exhibits the same tumbling statistics 
as a single dumbbell in shear flow. 

To compare the tumbling time r*"™ with the relaxation 
time r''°', we performed DPD simulations in a 18x18x18 
box and other parameters as before and followed the dy- 
namics of 6 dumbbells for 10000 time units. We computed 
their extension R{t) and azimuthal angle (j){t) every 3At. 
The intermittent nature of such time traces is immediately 
clear from the example shown in the inset of Fig. [3l 

In Table [D we list the average tumbling times r'"™ for 
different values of i?max and i?FENE- The values of r*"™ 
were computed from the angular dynamics as explained in 
the caption of Fig. [31 We have also determined the indi- 
vidual tumbling times from the time interval between two 
stretched configurations with a contracted one in between. 
Both methods yield comparable results, but contrary to 
single polymer experiments for which angular resolution is 
a limitation |26j . analysing the angular time traces works 
best numerically and is least dependent on the thresholds 
used to define a tumbling event. 

In our simulations, the distribution of tumbling times 
shown in Fig.[3|decays exponentially — this is known to be 
a robust property of a single polymer in shear flow [^^11^ . 
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T = time between two tumbling events 



Fig. 3: Distribution of tumbling times for 6 polymers in shear 
flow at Wi = 1 during 10000 simulation time units (provid- 
ing typically 2000 tumbling events). Dumbbell parameters are 
Hfene ~ 1 and Timax = 4. The distribution is well fitted by an 
exponential. Inset: Corresponding time trace of the azimuthal 
angle of a dumbbell. We take the tumbling time r to be the 
time between two large values of (j> (outside the dashed blue 
lines \(j>\ = 1.2) such that (j> takes at least one value close to the 
average (inside the grey region) for some intermediate time. 

The fact that there is only one intrinsic relaxation time in 
our model is another important illustration that its vis- 
coelastic behaviour is a good mesoscopic representation of 
the single-relaxation time in Oldroyd-B type constitutive 
equations. Qualitatively, the collective tumbling of DPD 
polymers is very much like what one expects from Brow- 
nian dynamics of a single dumbbell: the tumbling time 
increases with the fluid relaxation time (we almost have 
proportionality for i?max ==5)- As expected |22j . we find 
that r'"™ (X T™' for not too large Wi. 

We have also performed simulations with a very few 
dumbbells - only 12 out of the 8 • (18)'^ particles were 
paired up with FENE springs. Except for small string 
stiffness (iJpENE < 0.1) (presumably because of the poor 
statistics and the force singularity at large stretching) , cal- 
culation of the associated tumbling times gave values very 
close to those with 100% dumbbells listed in Tabled To- 
gether with the results of the previous section, this shows 
that the elasticity of the surrounding dumbbells docs not 
significantly affect local dynamics. It implies that by in- 
creasing the concentration of dumbbell polymers /p, we 
can increase the total normal stress without changing the 
properties of the fluid on a mesoscopic scale. 

Turbulent drag reduction. — To demonstrate that 
our approach can also be employed to study non-trivial 
macroscopic flows, we perform, for the first time to our 
knowledge, DPD simulations of high Reynolds number 
turbulent fiows in the presence of polymers. This sys- 
tem is known to exhibit drag-reduction - the phenomenon 
that the addition of minute amounts of polymer signifi- 
cantly suppresses turbulence and hence the drag [7]. It 
has recently been successfully captured by direct numeri- 
cal simulations [?fll^ and by studying the exact coherent 
states [ini- We have shown previously [TS| that the stan- 
dard DPD model reproduces the main characteristics of 
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Fig. 4: (a,b) Plots of the turbulent amplitude A{t) versus time 
as a function of Reynolds number Re for 10% (left) and 30% 
(right) polymer fractions. The Reynolds number Re is given 
besides each curve. The dotted line corresponds to the time 
at which we start switching off the driving force, (c) Influence 
of polymers on the stability of turbulent state in our dumbbell 
DPD simulations. As the fraction of polymers is increased, the 
onset of turbulence is shifted to higher values of Re and the 
turbulent amplitude is reduced. 

turbulent plane Couette flow (albeit the onset Reynolds 
numbers are somewhat too high perhaps due to the sig- 
nificant compressibility of the DPD fluid). We now study 
the effect of polymers on the turbulent state as a function 
of polymer fraction. 

To do so, we now switch to the soft FENE dumbbell 
interaction ^ for the reasons we explained before. The 
parameters for these runs are a ~ 4, p = A^T = 0.1,T = 1, 
^^FENE = 0.5 and Rc = 2 (the rheological properties are 
r]s ~ 0.69, rjp « 0.25/p and « 0.9); the simulation 
box has dimensions x Ly x ~ 60 x AO x 2cQ. 

We first drive the flow as in [15] with a force field that 
generates an array of streamwise vortices that are known 
to dominate the turbulent state at least close to the onset 
[50] . The driving field is turned on at time 200, then slowly 
ramped up to time 440, kept at a steady value till time 
1160, and then slowly turned off till time 1400. By follow- 
ing the maximum deviation A{t) from the linear flow field, 
A{t) = max(y(z) — jz) where V{z) is the spatially aver- 
aged mean streamwise velocity profile, it is then tested 
whether turbulence remains or whether it decays to the 
laminar state. The behaviour of A{t) was studied as a 
function of polymer fraction fp and the Reynolds number 
Re. We have varied the forcing amplitude over a factor 
12 and found the final turbulent state (if it exists) to be 
independent of the forcing amplitude. Typical data for 

^In these runs, the walls are modelled by a layer of continuous 
DPD matter (see 1151 ). rather than by a smooth potential. 



p-5 



E. Sultan et al. 




Y 40 Y 40 



Fig. 5: (colour online) Streamwise velocity and polymer ex- 
tension in the gradient-spanwise plane in a simulation with 
fp = 10% and Re = 2520. Colours represent the streamwise 
deviation (locally averaged) from (a) the laminar velocity (vec- 
tors denote the in-plane motion) and (b) the average dumbbell 
stretching. Colour code is such that red and blue indicate a 
±79% deviation from the average velocity, and for stretching 
a variation of ±66% about the average stretch (« 1.1). 

two polymer concentrations is shown in Figs. |3Jab). The 
saturated values of A{t^oo) are plotted in Fig. \M,c) as a 
function of Re for various fp. This stability diagram illus- 
trates the occurrence of drag reduction within our dumb- 
bell DPD model: for increasing polymer fraction, the onset 
of turbulence is pushed to higher Reynolds numbers and 
the turbulent amplitude is reduced. A similar stability di- 
agram was reported for the exact coherent states in plane 
Couette flow [29]. Finally, as Fig. [5] illustrates, not only 
do the dumbbells get preferentially stretched inside the 
streaks (blue and red regions in the velocity plot), as one 
might expect [15|, but the fluctuations in the stretching 
appear to be much larger than those in the velocity. 

Conclusion. — In this paper we have shown that the 
dumbbell DPD model has all the main characteristics of 
the Oldroyd-B class constitutive equations, while the sta- 
tistical properties of the dumbbells are similar to those of 
single polymers. It can thus be used to study polymer rhe- 
ology simultaneously at the macroscopic and mesoscopic 
scales. This makes the DPD dumbbell model an attractive 
candidate for simulations of polymer rheology in complex 
geometries and flows, in particular in situations where im- 
pact of the macroscopic flow on the mesoscopic statistical 
properties is of interest. We expect the model to be espe- 
cially useful and easy to implement for studies of cross-slot 
geometries, flows past objects, extrusion, etc. Detailed 
information about coarse grained polymer conformation 
readily available in our model may even shed new light on 
the mechanism of drag reduction. 
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